Adding a reaction-restoration type transmission rate dynamic-law to the basic SEIR COVID-19 model

The classical SEIR model, being an autonomous system of differential equations, has important limitations when representing a pandemic situation. Particularly, the geometric unimodal shape of the epidemic curve is not what is generally observed. This work introduces the βSEIR model, which adds to the classical SEIR model a differential law to model the variation in the transmission rate. It considers two opposite thrives generally found in a population: first, reaction to disease presence that may be linked to mitigation strategies, which tends to decrease transmission, and second, the urge to return to normal conditions that pulls to restore the initial value of the transmission rate. Our results open a wide spectrum of dynamic variabilities in the curve of new infected, which are justified by reaction and restoration thrives that affect disease transmission over time. Some of these dynamics have been observed in the existing COVID-19 disease data. In particular and to further exemplify the potential of the model proposed in this article, we show its capability of capturing the evolution of the number of new confirmed cases of Chile and Italy for several months after epidemic onset, while incorporating a reaction to disease presence with decreasing adherence to mitigation strategies, as well as a seasonal effect on the restoration of the initial transmissibility conditions.


Introduction
A novel Coronavirus (SARS-CoV-2) emerged from the city Wuhan in China in December 2019 and has caused a devastating public health impact across the world [1]. As of June 28, 2021, COVID-19 has caused over 180 million confirmed cases and over 3.5 million deaths worldwide [2]. The curves for daily confirmed new cases of COVID-19 in different countries present a high variability in their geometric forms. Every such curve shows a sequence of outbreaks and valleys when observed over time, while the sharpness of the outbreaks and the length of the valleys can vary [3].
In fact, in [3] we can find figures that show the daily number of confirmed new cases (7 day moving average), in particular from March, 2020, to January, 2021, for different countries. countries show a sharp first outbreak followed by a plateau of low height (that lasted several months) before the second wave, whose peak out-measures the peak of the first outbreak by at least three-fold. On the contrary, several South American countries presented an initial exponential phase of several months, soon after which a second wave of similar peak height as the first occurred. Finally, there are countries in which the curve of new cases of COVID-19 has shown extreme behaviors in some part of its evolution as compared to European and South American countries. For instance, Czech Republic has experienced a very weak first outbreak followed by a low plateau lasting for months; Iran did have a more pronounced first outbreak, but it was followed by a plateau of important height; and Indonesia experienced an almost unnoticeable first outbreak that is part of an exponential growth phase when zoomed out, which lasts several months. Classical compartmental models based on the classical Kermack & McKendrick SIR model [4] with constant parameters often used to model epidemics do not reflect the behavior over several months described above. Neither the (β, γ)SEIR model for a population of size Nwhich is compartmentalized into susceptible (S), exposed (E), infectious (I) and removed (R)given by S 0 ¼ À bSðI=NÞ; E 0 ¼ bSðI=NÞ À eE; I 0 ¼ eE À gI and R 0 ¼ gI; nor extensions of it have been efficient in adjusting the data well beyond the first epidemic outbreak when considering the transmission rate β and the removal rate γ constant. This is due to the unimodality of the active-infected-curve those models provide, i.e. one bell-shaped infection curve and an epidemic growth limited by the proportion of susceptible individuals [5,6]. In general, the epidemiological data series do not reflect that the percentage variation of susceptibles per proportion of active cases, i.e. |S 0 /S|/(I/N), is approximately constant, as is assumed in the aforementioned classical epidemiological models. In fact, there exists literature that evidences the changing temporal behavior of disease transmission in epidemic or pandemic situations [7][8][9][10][11][12][13][14][15]. In particular, there are studies using mathematical models-some aiming to understand COVID-19 transmission-that include the decrease in the transmission rate [7,9,[16][17][18][19][20], and some incorporating human behavioral factors as part of the cause for a temporal change in transmission. For instance, in [7] the behavior of the transmission rateproviding exponential saturation for a large number of infectives-for three consecutive months is shown for four geographical settings: worldwide, United States, Russia and Canada. For each setting, the trend is a decrease in the transmission rate for then stabilizing at a value several times lower than initially. They as well include a parameter representing mask wearing within their transmission rate. Finally, this study also shows that the recovery rate R 0 /I is for each setting much more stable than the transmission rate. Supporting the idea that the lack of efficiency for the classical compartmental models to adjust well to data is due to that the transmission rate β is assumed constant over time.
In this article we attempt to break the unimodality of the active-infected-curve of the classical epidemiological models. We introduce a novel way to model the behavior of the transmission rate β, considering a balance equation between a reaction rate and a restoration rate; and including the resulting dynamic law for the transmission rate into the classical SEIR model. The paper is structured in the following way: In the next section The Transmission Rate we provide some understanding about the transmission rate of infectious diseases. In the section The βSEIR Mathematical Model we introduce and describe a new basic model, which we call βSEIR model, by adding to the classical SEIR the aforementioned dynamic law for the transmission rate, and show some mathematical analysis. In the section Numerial Results we provide simulation results; and finally, the last section contains the discussion and conclusions.

The transmission rate
There exist at least two groups of epidemic control measures. The first, aims to reduce the population that is being hit by the disease, i.e., the susceptible population. Such measures are, for example, vaccination or limiting the mobility of individuals. The second, intents to reduce the force of infection that is defined as the product of three quantities: number of close contacts per unit of time of a susceptible individual (p C ), probability that a close contact is with an infectious individual (p I ) and probability of transmission given a close contact with an infectious (p T ). Reducing p I results in that per unit of time, there exist fewer active cases in the population, which is accomplished by eradication, i.e., removal from the system (e.g., slaughtering of sick animals which is widely used in animal epidemics, or banishing infectious individuals as was done aforetime), or by applying actions for a rapid recovery. Notice that the product of p C and p T is called the transmission rate and is usually denoted by β (see e.g. Eq (1)). Hence, the objective of most mitigation strategies that aim to reduce the force of infection, aim to reduce the transmission rate (lower β) by either increasing physical distance and hence reducing the number of close contacts (lower p C ) or blocking the transfer of pathogens to a new host (lower p T ). There are secondary measures such as for example reducing population movement (which is not reducing physical distancing nor blocking transmission), which make close encounters less likely.
When a highly transmissible disease with high mortality or morbidity invades a population of mostly susceptible individuals, and a vaccine is not in sight in the short term (as was initially the case for COVID-19), health authorities' only way for reducing morbidity and mortality is mitigation, while the general population's duty is to comply to the new norms and desired behavior. In other words, the efforts are put into reducing the transmission rate β. In this sense, β is a time-spacial dependent variable, i.e. it changes according to time and location. Also less controllable physical-environmental aspects in relation with the bio-chemical characteristics of the pathogen may influence it (but we consider those factors constants in this study). It is also worth mentioning that in general populations, individuals may live and participate in several cultural regions, which may also determine the variability of the transmission rate. We will assume in this study that the population stays within its territory during the time horizon studied, behaving homogeneously in this respect. We suppose that the disease studied will have a base line transmission rate β 0 , that we will call natural transmission rate, measured for a population that is initially free from the disease and does not consider any mitigation strategies or personal protective measures.
One of the characteristics of COVID-19 was that it has had a large media coverage since the first confirmed cases appeared in December 2019. This provoked sentiments of fear in the general population and played an educational role for pandemic preparedness (e.g. global media emphasizing on washing hands and physical distancing) before the imminent arrival of COVID-19 in many countries. It shaped how countries would confront COVID-19 right from the appearance of their first confirmed case and even before. The reproduction numbers corresponding to different geographical locations was most likely to be between 2 and 5 [21,22], and was shown to rapidly decrease during the first weeks of the pandemic [12][13][14], however reaching values above one that nevertheless allowed COVID-19 expansion.
To consider this decreasing effect, many authors assume an exponential decay of the transmission rate for a certain amount of time, for example, in [9,18] they assumed in their continuous time model β(t) = β 0 exp(−b 0 t), t � 0, or in in [19,20] they defined β k = β 0 a k , 0 < a < 1, k � 0, where k is a day-counting integer in their discrete time model. In order to extend the horizon of validity of their model some authors consider an exponential decrease from β 0 to a minimum positive bound [15,23]. To understand the rate of decrease from that baseline natural β 0 and identify a time varying β(t) transmission rate, some researchers use mathematical expressions and the data on active cases I(t) and removed cases R(t) in a population of constant size N; for instance, one can use β(t) = −NS 0 /(SI) as in [7], which can be obtained from a SIR model and approximating the derivatives using the finite differences on one week running averages; or one could use β(t) = γ + I 0 /I at the beginning of an outbreak, assuming S * N in the SIR model.
More time-varying transmission rates have been considered within mathematical models. For instance, the authors in [9] capture the early decreasing trend of COVID-19 in Malaysia using a time varying exponential decay log function β(t) = zβ(1 − p) t for the transmission rate in an SIR model, that uses a fractional term z to measure the effectiveness of interventions and a proportion p to account for depletion. In the literature there are studies that, in order to capture realistic disease transmission, assume non-linear functions of S and I governing the force of transmission, as for instance in [16], where the force of transmission they use in an SIR-type model depends on the product of fractional powers of S and I. They use the model to fit COVID-19 data of Italy, Germany, France and Spain.
In [24] the authors include in an SIR model time-varying transmission rate, assuming that the probability of transmission of a susceptible is βλ t (I t /N), where λ t (�) is a random variable, they refer to this model as a Spatial-SIR model. In [25] the authors assume a contagion rate as a sum of a base-line transmission rate and a component that satisfies a first order linear differential equation to represent the effect of non-pharmaceutical interventions (NPIs).
Behavioral factors represented in the transmission rate are also considered by more authors in order to represent the changing dynamics of the transmission rate. For instance in [26] the authors study a model that incorporates a non-constant transmission rate β(M) that depends not only on the current number of infectious individuals but on M, representing an information index that summarizes the current and past history of disease prevalence. Part of their results discuss that social behavioral change may trigger oscillations. The study in [27] extends an SIR type model defining a transmission rate that captures the impact of school and workplace closure through a function of time. The changing behavior of this function is based on Imitation Dynamics [28] and describes population-level support dynamics for closure. The article in [29] also uses Imitation Dynamics and studies a population in which individuals develop and learn a behavior of mutual protection.
The novelty of this article is that, at each moment in time we consider the variation of the transmission rate to be given according to a balance equation between two opposite thrives: a reaction rate and a restoration rate.
In what follows we are going to justify the functional forms of the reaction and restoration rates, as well as present the βSEIR model that incorporates the dynamic law for β. Further, we are going to analyze its effect on the shape of the main epidemiological curves.

The βSEIR mathematical model
In this section we present the dynamic law of the transmission rate β in order to introduce the βSEIR model and some mathematical analysis.

Transmission rate β
In this subsection we derive the form of the transmission rate at which susceptible individuals become exposed upon contacts with infectious. The only infectious class of the model is the I class. We model the case of an infectious disease transmitted directly from person to person, and assume that at the beginning of an outbreak, the appearance of first cases do not provide a reason for alarm and panic. Hence, initially the disease propagation is due to a high natural-transmission rate intrinsic to the population while no interventions to mitigate disease spread are in place. We call this natural-transmission rate β 0 . In general, β 0 makes the disease expand rapidly, producing a fast initial increase in new cases.
We present a new form for the transmission rate, which is represented by a dynamic, timedependent quantity that is governed by a balance equation between a reaction rate, g[t, I], and a restoration rate, f[t, β], i.e. the proportion by which the transmission rate decreases and increases per day, respectively, represented by Next we justify the introduction of a reaction rate and a restoration rate and propose a functional form for each.
Reaction rate. During severe epidemic outbreaks that attract huge public attention and media coverage due to for instance a high morbidity and/or mortality in the population, a steady increase in the implementation of measures that aim to reduce the transmission rate can be observed. As long as there is no licensed vaccine or treatment, these measures are mainly based on non-pharmaceutical interventions. Here we are interested in those directed to diminish the factors p C and p T , whose product defines β, such as social distancing measures (that reduce the number of close contacts between people: large-scale or home quarantines, workplace non-attendance, travel restrictions, prohibition of social gatherings, school closures, etc.) that reduce p C , or blocking measures (that, given e contact, reduce the pass of the pathogen: hand-washing, respiratory etiquette, face-masks usage, etc.) that reduce p T , see [30,31]. When fear governs the population, people react, complying with mandatory measures or adopting self-protecting measures to avoid infection [32]; we note that risk communication is also a factor to support the general public response [33]. It is to expect that the higher the severity of the disease is, the more effort is put into mitigation. In this sense, individuals' reactions produce a decrease in the transmission rate from its initial natural-transmission rate value, β 0 . Hence, we define a reaction rate that we denote g[t, I], which is non-negative and positively correlated with the number of active cases I(t), in a way that it increases when I(t) does. It may also depend on other circumstantial conditions of the moment relative to the population. These we will discuss further later in the text.
We assume in this article, that the reaction rate, follows the Michaelis-Menten model [34] describing the reaction to the presence of the infectious (active cases) I(t) at any moment in time, i.e. we define where we call μ(t) the reaction coefficient, which is a non-negative function that represents the daily maximum possible reduction at time t, and I m > 0 is the half-saturation constant, i.e. is the number of active cases, where the reduction is half-maximal. Notice that the parameter I m characterizes the population, i.e. it determines its sensibility to react to active cases. Restoration rate. It is important to observe that upon the appearance of a reaction rate there exist socio-environmental factors that tend to restore the transmission rate to the level observed at the beginning of the pandemic, i.e. to β 0 [27,35]. When in a certain location the health authorities cease to impose protective measures, e.g. the use of face masks, and individuals lost their initial fear, then the transmission rate that had been reduced due to these measures no longer stays low and returns to its natural level. Therefore, we introduce a restoration rate that at each instant t, t > 0, is responsible for an increase in the transmission rate. The restoration rate that we denote by f[t, β], is a non-negative function that correlates directly with the distance between β(t) and its natural value β 0 , as presented in the following equation: where a 2 R is a positive exponent. We call ν(t) the restoration coefficient, which is a non-negative function that regulates the daily form of the restoration rate. Also, f is an increasing function of |β(t) − β 0 |, and f[t, β 0 ] = 0, which means that if the transmission rate β(t) reaches at a certain time point t its natural value, then there is no deviation to restore.

The βSEIR model: Formulation and analysis
We incorporate the differential Eq (2) to the classical SEIR model (1) obtaining the βSEIR model, given by the following system of equations with some non-negative initial conditions I] as in Eqs (4) and (3) respectively. Table 1 describes the variables and parameters of the model. In the following we show that 0 � β(t) � β 0 , for all t � 0 and β 0 > 0. Observe that β 0 � − μ (t)β, and hence due to Grönwall's inequality (see [36]) we can conclude that bðtÞ � b 0 e À R t 0 mðsÞ ds . Therefore, β(t) is non-negative for any non-negative initial condition β 0 .
We also observe that β(t) � β 0 is an equilibrium solution of the first equation in system (5) as long as I(t) � 0. In case I(0)>0, i.e. disease is present in the population, β 0 (0) < 0 and hence β

PLOS ONE
(t) decreases initially. Since β 0 � ν(t){|β − β 0 |/β 0 } α β, using Grönwall's inequality and the fact that the solutions to the differential equation β 0 = ν(t){|β − β 0 |/β 0 } α β that pass through points ðt;bÞ with 0 <b < b 0 are increasing and bounded by β 0 , we obtain that β(t)�β 0 for all t � 0 and β 0 > 0, I 0 > 0. Just as for the classical SEIR model, we observe that the epidemiological state variables of the model in system (5) remain positive for positive initial conditions and are bounded by the total population size N. Also, adding the second, third and fourth equations in system (5) together we obtain (S + E + I) 0 = −γI < 0 whenever I > 0. Hence, S + E + I is a non-negative smooth decreasing function, and therefore lim t!1 ðS þ E þ IÞ exists. On the other hand, the derivative of any smooth non-negative decreasing function must tend to zero, and hence 0 ¼ lim Similarly, by adding the second and third equation in system (5) together, we can prove that lim exists when time tends to infinity, we can then conclude that lim Hence, the long term behavior of the model we present holds (i.e. the limit of the epidemiological state variables exist for infinite time), just as for the classical compartmental epidemic models SIR or SEIR with constant transmission rate β [37]; in particular we have shown that the disease in the long term goes extinct. Finally, as for classical models we can obtain a threshold condition that determines an initial epidemic outbreak. Given β 0 > 0 and considering that β(t)�β 0 for all t � 0 as well as that the state variables are bounded by N, we have that ( Hence, we define the basic reproduction number [38,39] as R 0 ≔ b 0 =g. This way, if R 0 < 1, then the curve representing the infectious population is decreasing to zero and there is no epidemic outbreak. On the other hand, if R 0 ¼ b 0 =g > 1, then (E + I) increases initially when we assume S * N and β(0) = β 0 , and increases as long as (β(t)/γ)S(t)/N > 1 holds. Notice that, in this case, the curve of infectious individuals may be increasing at several time intervals depending on the behavior of the function β(t), and not only depending on the ratio of susceptible individuals in the population that is decreasing according to the second equation in system (5). In Section 1 we will call the quantity (β(t)/γ)S(t)/N the Effective Reproduction Number, R e , and its dynamics will determine disease dynamics.
In the context of our study, we consider short-medium term scenarios under an epidemic situation, i.e. when R 0 > 1. We also consider no replacement of susceptibles, due to the time frames we describe and assuming that: first, infected individuals acquire some kind of protective immunity for several months after infection [40]; second, we do not consider multiple variants that may provoke new infection. In particular, we will point out the differences compared to the classical SEIR model with constant β, in which E and I are variables that describe the known unimodal behavior of one bell-shaped curve.
Our model generalizes an idea presented in [41], where the authors consider an SIR-type model with variable transmission rate of the form β(D(�)), in which the time dependent function D(�) represents social distancing that individuals in the population maintain to each other, governed by the dynamics of D 0 = −λ 1 (D − D�) + λ 2 I/N, where λ 1 and λ 2 are positive constants, and D � is the culturally dependent natural social distance, to which D(�) converges in the absence of disease (I = 0).

Numerical results
In this section, we present simulations of disease dynamics assuming: first, constant reaction (μ(t)) and restoration (ν(t)) coefficients; second, we extend our results to consider time-varying reaction coefficients, representing a diminishing mitigation effect of government responses; third, an additional seasonal effect on the restoration coefficient and we show that our model is capable of capturing real COVID-19 disease data. For the simulations we use the Python programming language [42].

The autonomous βSEIR model: Constant reaction and restoration coefficients
We consider in this subsection the autonomous βSEIR model from system (5), assuming a constant reaction coefficient, μ(t) = μ; i.e. we suppose that the reaction to reduce the transmission rate just depends on point prevalence levels and not on an additional time factor (see Eq (3)); and also assuming a constant restoration coefficient, ν(t) = ν; i.e. the regulation on the restoration rate depends only on the deviation of the transmission rate from its natural value β 0 and not on external temporal factors (see Eq (4)). We present qualitative results of our model through simulations considering an 18 months time span, and present in each of the Figs 1-3 five subplots that represent the dynamics for: restoration f(t, β) and reaction rates g(t, I); the transmission rate β(t); the effective reproduction number R e ; the new confirmed cases eE(t); and finally we illustrate the cumulative number of cases. The effective reproduction number is a time-varying threshold quantity-defined by R e ðtÞ ≔ ðbðtÞ=gÞSðtÞ=N-such that the number  Table 2.  Table 2.

PLOS ONE
of cases increase while R e ðtÞ > 1, reach a peak when R e ðtÞ ¼ 1 and decrease when R e ðtÞ < 1 [43,44], and in particular R e ð0Þ ¼ R 0 . Additionally, the constant restoration coefficient ν takes in Figs 1-3 the values 0.8 (high), 0.5 (medium) and 0.2 (low), respectively, representing high, medium and low rates to restore transmission levels due to the urge to return to the natural transmission rate β 0 . For each constant restoration coefficient value, we choose within each figure the constant reaction coefficient μ to take the values 0.3 (low, in blue), 0.4 (medium-low, in green), 0.5 (medium-high, in orange) and 0.6 (high, in red), representing different constant levels of the daily maximum reaction to disease that reduces the transmission rate. The remaining parameter values and initial conditions used in the figures are described in Table 2. Note that the basic reproduction number obtained from the values β 0 = 0.65 and γ = 1/14 from Table 2 is R 0 ¼ 9:1 > 1. For illustration purposes, we use a value significantly larger than one. Fig 1 considers a high restoration coefficient, ν = 0.8. Observe that for high reaction coefficients, e.g. μ = 0.6 as well as for small reaction coefficients, e.g. μ = 0.3, the restoration rate f(t, β) and the reaction rate g(t, I) are very similar to each other. Initially the restoration rate is slightly smaller than the reaction rate, producing that β 0 (t)<0 (from the first equation in system (5)) and therefore being the reason for the initial decrease of in the transmission rate. After that drop, the reaction and restoration rates are almost equal, producing a plateau in the transmission rate due to β 0 (t)*0, and subsequently, a slightly larger restoration rate than reaction rate produces an increase in the transmission rate, which converges to its original natural transmission rate value β 0 . Additionally, observe that for all values of the reaction coefficient μ, the duration of the plateau in the transmission rate is directly correlated with the value of μ: the higher the μ value, the longer its duration. Note that, while the transmission rate is at the plateau-i.e. β(�) behaves similar to constant-when equaling the first equation in our βSEIR model (system (5)) to zero, we conclude from f(t, β) = g(t, I) when α = 1 that ; with l ≔ ðn À mÞ=n: Hence, if I m � I, then β * β 0 λ. Indeed, notice that the value of β 0 λ for μ = 0.3, 0.4, 0.5 and 0.6 in Fig 1 are respectively 0.406, 0.325, 0.244 and 0.163, which correspond very closely to the plateau levels of the respective transmission rates. For all μ cases, the effective reproduction number shows a decreasing shape, staying above one at least during the first months of a pandemic (the horizontal dotted line represents R e = 1 in the figure). While the effective

PLOS ONE
A reaction-restoration type transmission rate dynamic-law within the basic SEIR model reproduction number stays above one, one can observe that the number of new confirmed cases increases, reaching a peak when the effective reproduction number reaches the threshold value one. The lower the reaction coefficient value is, the larger is the transmission rate throughout the epidemic, which produces sooner and larger epidemic peaks of new confirmed cases, as well as a rapid increase in the number of cumulative cases, reaching quickly a number close to the final number of infected individuals throughout the whole epidemic. This happens shortly after the effective reproduction number reaches the value of one, and is due to a small number of susceptibles remaining in the population at that time. Fig 2 considers a medium restoration coefficient, ν = 0.5. We observe that for a high reaction coefficient μ = 0.6, the restoration (f(t, β)) and reaction (g(t, I)) rates show an initial oscillation and differ from each other more clearly, being initially the reaction rate larger than the restoration rate and subsequently both intersecting several times. This produces an oscillatory behavior in the transmission rate due to the β equation in system (5), attaining the transmission rate a local maximum or minimum value each time the reaction and restoration rates intersect, i.e. f(t, β) = g(t, I). On the contrary, in the case of a low reaction coefficient value, e.g. μ = 0.3, we observe very similar restoration and reaction rates (blue curves) just as in Fig 1. We can also observe from Fig 2 that a high (red) or medium-high (orange) reaction coefficient μ, drives the effective reproduction number below one much faster than in was observed in Fig 1 (in the case of a higher restoration coefficient) and way before the cases for medium-low (green) and low (blue) reaction coefficient values. It is interesting to see that in the cases of higher reaction coefficients (red and orange), after the initial fast drop of the effective reproduction number, during the remaining time pictured it oscillates around one. When comparing the curves of the new confirmed cases and cumulative cases for these two reaction coefficient values, with the cases of low and medium-low reaction coefficients, we notice that the number of cases is way higher for the latter. Additionally, one bell-shaped curve of new confirmed cases is being observed for smaller reaction coefficient values (blue, green), since the effective reproduction number only manages to cross the threshold R e ¼ 1 (see dotted line) once due to the small number of susceptibles remaining after that first large peak. On the contrary, an oscillatory behavior is seen for higher μ values, obtaining small epidemic peaks each time the effective reproduction number reaches one in a decreasing manner.
In other words, the unimodality of the behavior for new confirmed cases observed when μ is low (blue) or medium-low (green) is broken for high (red) and medium-high (orange) μ values, which was not seen in Fig 1. One can also observe that for higher μ values (red and orange), the increase in the cumulative cases is close to linear in the time-frame pictured, as opposed to the rapid increase in cumulative cases for smaller reaction coefficient values (blue, green). Fig 3 shows the case of a low restoration coefficient ν = 0.2. We observe oscillatory behavior in the restoration and reaction rates, producing an oscillatory behavior in the transmission rate, and hence also in the effective reproduction number and in the number of new confirmed cases. We observe that for high reaction coefficient values, e.g. μ = 0.6, the absolute difference between the reaction and restoration rates are larger and their intersections occur sooner, and hence the oscillations in the transmission rate have a higher amplitude and their local maxima occur sooner, than in the case of lower reaction coefficient values, e.g. μ = 0.3. Also, transmission rates with higher amplitudes produce less pronounced peaks in the oscillatory behavior of the number of new confirmed cases. Additionally, sooner starting oscillations correspond to higher values of the reaction coefficient μ. The cumulative cases show a near to linear increase, where smaller slopes correspond to higher reaction coefficient values.

A non-autonomous βSEIR model: Time-varying reaction coefficient representing a diminishing mitigation effect
It is to expect that, when a disease enters a population, the reaction coefficient μ(t) right after the onset increases quickly, reaches a maximum value μ 0 , and then decreases due to many factors. This can be deduced (at least) from two sources: (a) The data and information provided by the Oxford COVID-19 Government Response Tracker (OxCGRT) and the time curves defined by the Stringency Index [45,46]. This index records the intensity of several government responses combined, such as containment and closure policies by country. In general, we observe that, first, the time curves representing the stringency index rise, but then follow a decreasing behavior due to local socioeconomic reasons [47][48][49][50][51][52]. (b) Studies in behavioral science explain the public fall of adherence to mitigation (distancing) measures, as a daily average compliance curve in times of COVID-19 shown for instance in [53] or discussed in the conclusions in [54]. Additionally, it is known that information-based interventions positively impact compliance with mitigation restrictions, such as keeping a certain social distance, decreasing the number of times individuals go out and their time spent outside. Nevertheless, if individuals have been restricted for a prolonged period of time, compliance with such mitigation strategies decreases [55]. In this manuscript we refer to a time period right after disease onset, where there is no vaccination yet available.
We can include such a behavior-representing a diminishing mitigation effect of restrictive measures-through a time-varying reaction coefficient μ(t), for instance, given by the following Eq (6), with 0 < μ 0 < 1. Notice that h(�) is a non-negative, unimodal function such that h 0 (b) = 0 and h(b) = 1, i.e. it achieves its maximum value at t = b and then decreases at a rate that depends on the parameter a.
In each of the Figs 4-6 we show for high, medium and low constant restoration coefficients ν, respectively, the curves for the restoration (f(t, β)) and reaction (g(t, I)) rates, the transmission rate β(t), the effective reproduction number, the new confirmed cases and the cumulative cases, for a forgetting curve h(t) of the form given in Eq 6, with a = 40, b = 90, such that the maximum occurs at t = 90 days. Within each plot we present curves for different maximum values μ 0 of the now variable reaction coefficient μ(t) = μ 0 h(t). The other parameters used for these figures are given in Table 2.
We present in Fig 4 the case of a high constant restoration coefficient ν = 0.8. We observe that the restoration and reaction rates are very similar regardless of the μ 0 value. Despite their similarity, initially, up to day 90, a slightly higher reaction than restoration rate produces a sharp decrease in the transmission rate. The decreasing shape of the forgetting curve starting on day 90, immediately reduces the reaction rate to slightly below the restoration rate, producing an increase in the transmission rate, eliminating the plateau observed in Fig 1, in which the reaction coefficient was assumed constant. Due to the incorporation of a forgetting-curve in the reaction coefficient we can also observe higher and earlier occurring epidemic peaks in the curves of new confirmed cases as compared to the constant reaction coefficient case (see Fig 1), especially for large μ 0 values. Fig 5 depicts the dynamics for a medium restoration coefficient ν = 0.5. Here we can observe-especially for a large maximum value μ 0 of the reaction coefficient μ(t) = μ 0 h(t)-that the restoration rate f(t, β) and the reaction rate g(t, I) differ more from each other, as compared Simulations for a high restoration coefficient ν = 0.8. The first subplot depicts the shape of the compliance curve h(t) = (a + 2b)t/(t 2 + at + b 2 ), with b = 90, a = 40, and the second plot in the first row illustrates the restoration rate f(t, β) (dotted) and the reaction rate g(t, I) (solid). The remaining subplots show: The transmission rate β(t), the effective reproduction number R e ðtÞ, the new confirmed cases eE(t), and the cumulative cases E(t) + I(t) + R(t). The maximum value μ 0 of the reaction coefficient (μ(t) = μ 0 h(t)) used in each subplot is: μ 0 : 0.3 (blue); 0.4 (green); 0.5 (orange) and 0.6 (red), and the remaining parameter values and initial conditions are as in Table 2.

PLOS ONE
to the case depicted in Fig 4. If we observe the curves for μ 0 = 0.6 we see that initially, up to day 90, the reaction rate is clearly higher than the restoration rate, which again produces a sharp decrease in the transmission rate. Every time the restoration and reaction rate curves intersect, we can observe a local maximum/minimum in the transmission rate, which eventually increases and converges back to its natural value β 0 . The changing human behavior reflected in the restoration and reaction rates, observed for μ 0 = 0.6, produces a late occurring but large epidemic peak, preceded by one small peak. This dynamics can also be explained by observing the shape of the effective reproduction number, which crosses the threshold R e ¼ 1 three times (red curve). Fig 6 shows the dynamics for a low constant restoration coefficient ν = 0.2. Here we observe initial oscillatory dynamics for the restoration and reaction rates, which explain the oscillatory dynamics in the transmission rate and also in the effective reproduction number, which crosses the threshold R e ¼ 1 (see dotted line) several times, producing small oscillations in the new confirmed cases, followed by a large epidemic peak.

A non-autonomous βSEIR model and the example of COVID-19 in Chile and Italy: Time-varying reaction and restoration coefficients
Although it has not been proven that the virus SARS-CoV-2 is seasonal in nature, seasonality of transmission may be an important factor to consider since social behavior is environmentally driven [56]. The main hypothesis regarding seasonality indicates that the higher the temperature the fewer infections occur [57,58].
Therefore, additional to a time-varying reaction rate, the restoration rate could depend, among other, on environmental factors such as for instance temperature. For example, in  Table 2.
https://doi.org/10.1371/journal.pone.0269843.g006 regions with predominantly Mediterranean climate, an annual oscillation of the daily mean temperature can be observed [59], as is the case in a large part of Italy or central Chile. The seasonality could also represent peoples' behavior during winter months vs summer months, or during vacation periods and special holidays. These factors affect the tendency of the population to return to its natural transmission rate β 0 , while making it seasonal in nature.
Hence, to capture seasonal factors we consider an annual periodic restoration coefficient of the form where ν 0 is the average restoration coefficient value, σ the amplitude of the oscillation and t max the moment in time, where the restoration coefficient is maximum. In general, social studies would be needed to determine the nature of the seasonality depending on peoples' behavior and habits in specific cultural and geographical contexts, independent of disease presence. This could help determine the parameters present in the seasonal part of the restoration coefficient. For instance, t max could happen during winter (where individuals meet inside), with a larger amplitude σ if a holiday coincides with cold weather and thus people tend to gather and restore normal conditions. Also, σ * 0 implies ν(t)*ν 0 constant, and hence there would be no seasonal variation in the restoration coefficient.
Considering this non-constant seasonal behavior for the restoration coefficient ν(t) as given in Eq (7), as well as a time-varying reaction coefficient as described in Eq (6), we illustrate through Figs 7 and 8 the applicability of our model during the first months of the pandemic, fitting the model to COVID-19 data of new confirmed cases of Chile and Italy [60,61]-minimizing in the least square sense the residuals between the outcome of our model and the data set-respectively from March 16th, 2020, to February 16th, 2021, and from February 24th, 2020, to October 31st, 2020. Table 3 shows the respective fixed and fitted parameter values used in each figure. Chile and Italy served as examples of countries located at different hemispheres that experienced COVID-19 in distinct ways, due to cultural differences, distinct levels of initial knowledge of the virus, seasons, etc., and this way permitted us to show that our model can capture different COVID-19 dynamics and partially describe them through reaction and restoration thrives of the population.
We can observe how the restoration rate and reaction rate differ among countries: In the case of Chile (see Fig 7), they are very similar, expecting a population whose reaction to disease Simulation fitting COVID-19 data from Chile. The first subplot depicts the restoration rate f(t, β) (dotted) and reaction rate g(t, I) (solid); the second plot the transmission rate β(t); and in the third plot the blue dots represent data of daily confirmed new cases of COVID-19 in Chile, from March 16th, 2020, to February 16th, 2021 [60]. The red curve represents the least square fit of the model to the data with parameter values as in Table 3 for the population of Chile, with N = 18 million individuals and initial conditions of the model (5)

PLOS ONE
presence and urge to return to their normal behavior (restoration) is similar in nature. Nevertheless, a small difference in those rates produce a large impact on the transmission rate, such as the steep decrease observed initially.
On the other hand, from Fig 8 one can see that the reaction rate (solid curve) and the restoration rate (dotted curve) differ more than in the case of Chile, having initially a population, which is reacting to disease presence faster than their urge to restore normal conditions, for later reversing their behavior three times, producing during that time period a long lasting plateau. Finally, one can see that a larger restoration rate than reaction rate produces the appearance of a large second peak.

Discussion and conclusions
In references from 1973 [69] and 1989 [70], strategic models are defined as those that, despite not having high resolution concerning a specific reality, have the advantage of containing all the minimum aspects of the referenced system. The main contribution of this work is to achieve a dynamically richer low-cost model, that is, one that adds only one more differential law to the classical SEIR model (without introducing new compartments in the population). Simulation fitting COVID-19 data from Italy. The first subplot depicts the restoration rate f(t, β) (dotted) and reaction rate g(t, I) (solid); the second plot the transmission rate β(t); and in the third plot the blue dots represent data of daily confirmed new cases of COVID-19 in Italy, from February 24th, 2020, to October 31st, 2020 [61]. The red curve represents the least square fit of the model to the data with parameter values as in Table 3 for the population of Italy, with N = 60.5 million individuals and initial conditions of the model (5)

PLOS ONE
There is evidence in the literature that populations change their behavior when facing dangerous diseases, i.e., reacting to these by managing to modify the transmission rate. Our proposed model, the βSEIR model, is a serious candidate to contain the minimum aspects for disease transmission of a high impact infectious-contagious disease in populations that, while living with the urge to restore normal conditions, react to reduce favorable conditions for disease transmission. Additionally, the latter occurs in a setting where vaccines are not available (and hence the model just considers one susceptible class), the disease is new to the population, and there are no multiple variants. Specifically, the novelty of the βSEIR model we present is, that it incorporates a variation in the transmission rate, which occurs proportional to the transmission rate itself but also proportional to a tension given by the difference (f − g) between a) reaction rate (g) to disease presence that may include other behavioral factors, such as compliance with mitigation strategies, and b) restoration rate (f) that aims to restore a certain intrinsic value of disease transmission, due to for instance socio-environmental elements.
Our results show an important gain in dynamic possibilities even in the case where Eq (2) in the βSEIR model given in Eq (5) is autonomous, i.e. f and g do not depend explicitly on time. Indeed, we can see in Figs 1-3 the appearance of several epidemic peaks and initial oscillatory dynamics, explained by the tension between reaction and restoration thrives of a population. In particular, we observe that high restoration coefficients ν-affecting the restoration rate f and representing a higher urge of the population to return to normal conditions-induce temporary stabilization of the transmission rate after an initial drop, being the duration of this plateau larger, the larger the reaction coefficients μ (affecting the reaction rate g) are, i.e. the higher the self-protective reaction is to disease presence. When considering small ν values (a small urge to return to normality), we observe oscillations in the transmission rate, with higher amplitudes for higher μ values, i.e. amplitudes are higher if individuals' reaction to disease presence is higher. These oscillations in the transmission rate generate oscillations in the effective reproduction number, which lay around one for a period of time proportional to the value of the reaction coefficient μ. Regarding the curve of new cases, we show that in general higher restoration coefficients ν produce unimodal behavior, whereas lower ν values generate the appearance of a finite number of peaks with decreasing peak size, in a way that for large reaction coefficients μ the timing between peaks is smaller. These results already define interesting future mathematical challenges.
In general, our results show how the transmission rate is impacted by the reaction rate g and the restoration rate f. In particular, we observe that a small difference between reaction and restoration rates may produce a large impact in the transmission rate. This highlights the importance of individual behavior in a pandemic setting, where even the behavior of a small number of individuals could change the dynamics of a disease drastically.
Curves of new confirmed cases with two and up to three waves in a one year time frame have been observed, differing among countries in time and size, with peaks and valleys of different heights, proper to a pandemic still under development in a population without protective immunity [3]. We notice that the βSEIR model can capture such patterns at the cost of varying the reaction coefficient μ(t) and for further dynamic richness varying the restoration coefficient ν(t). We can justify a time-varying reaction coefficient μ(t) that considers two aspects: first, it follows the shape of the Stringency-Index [45]-that records for each country government mitigation measures-, second, it reflects the reduced compliance with or adherence to mitigation strategies observed with time. The βSEIR model shows even richer dynamics when introducing such a time-varying reaction coefficient (see . Additionally, our model is capable of capturing the time series of new confirmed cases of Chile and Italy when including-additionally to including a time-varying reaction coefficienta time-dependent seasonal variation in the restoration coefficient ν(t), reflecting distinct temporal and possibly behavioral characteristics of two countries located at different hemispheres during their first pandemic year. Our results are, at first glance, good indicators for the richness that a model of such low structural complexity, as the one proposed here, can provide.